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This study demonstrates that coupling of a material thermal response code and a flow solver with non- 
equilibrium gas/surface interaction for simulation of charring carbon ablators can be performed using an 
implicit approach. The material thermal response code used in this study is the three-dimensional version 
of Fully Implicit Ablation and Thermal response program, which predicts charring material thermal 
response and shape change on hypersonic space vehicles. The flow code solves the reacting Navier-Stokes 
equations using Data Parallel Line Relaxation method. Coupling between the material response and flow 
codes is performed by solving the surface mass balance in flow solver and the surface energy balance in 
material response code. Thus, the material surface recession is predicted in flow code, and the surface 
temperature and pyrolysis gas injection rate are computed in material response code. It is demonstrated 
that the time-lagged explicit approach is sufficient for simulations at low surface heating conditions, in 
which the surface ablation rate is not a strong function of the surface temperature. At elevated surface 
heating conditions, the implicit approach has to be taken, because the carbon ablation rate becomes a 
stiff function of the surface temperature, and thus the explicit approach appears to be inappropriate 
resulting in severe numerical oscillations of predicted surface temperature. Implicit coupling for 
simulation of arc-jet models is performed, and the predictions are compared with measured data. 
Implicit coupling for trajectory based simulation of Stardust fore-body heat shield is also conducted. The 
predicted stagnation point total recession is compared with that predicted using the chemical equilibrium 
surface assumption. 


Nomenclature 

= mass fraction 
= diffusion coefficient, m 1 2 /s 
h = enthalpy, J/kg 

k t = thermal conductivity of translational temperature, W/m-K 

k v = thermal conductivity of vibrational temperature, W/m-K 

m = mass flux, kg/m 2 -s 

M = molecular weight, kg/kmole 

p = pressure, Pa 

q conv = convective heat flux at surface, W/m 2 
Vcond = In-depth conductive heat flux at surface, W/m 2 

q m = radiative heat flux at surface, W/m 2 

r c = corner radius, cm 

r n = nose radius, cm 

t = time, s 

t* = time of the last updated interface boundary conditions, s 

T = temperature, K 

T t = translational temperature, K 

T v = vibrational temperature, K 

T cf0 = environment temperature, K 

v = mass injection velocity, m/s 
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= source term of gas-surface interactions, kmole/m 2 -s 
= absorptance 

= maximum allowed variation on the surface temperature, K 
= emissivity 
= total density, kg/m 3 
= Stefan-Boltzmann constant, W/m 2 -K 4 

= time interval of updating the interface boundary conditions, s 
= char 

= pyrolysis gas 
= gas species 
= wall 


I. Introduction 

Vehicles designed for Earth entry at super-orbital velocities, as well as those designed for ballistic entry at 
orbital velocities, typically use thermal protection system (TPS) materials that pyrolyze and ablate at high 
temperature for mass-efficient rejection of the aero thermal heat load. For design and sizing of ablating TPS 
materials, it is imperative to have reliable numerical procedures which can accurately compute the vehicle 
aero thermal environment, the material surface ablation, and the material internal thermal response. It has been 
demonstrated that accurate prediction of ablative heat flux also requires a fluid-solid shape change coupling 
simulation. 1,2 In our previous fluid-solid shape change coupling computations, the ablative carbon surface was 
assumed to be at chemical equilibrium. Chemical equilibrium is a good assumption for many space entry 
applications, but it may not be valid for all conditions. Besides, the aerothermal environments computed by Data 
Parallel Line Relaxation method (DPLR) 3 were those for a fully catalytic non-ablating surface. Thus, an engineering 
correlation with blowing reduction parameter had to be introduced in the Two-dimensional Implicit Thermal- 
response and AblatioN program (TITAN) 4 simulation to take into account the effect of mass injection on reduction 
of convective heat flux. A mass transfer coefficient also had to be defined based on the heat transfer coefficient 
under the assumption of a relatively thin boundary layer for the computation of char recession rate. In our recent 
study 5 , integrated fluid-material response analyses with finite rate gas-surface interactions were performed. The 
coupling procedure, considering surface thermal chemistry and shape change, was achieved using an explicit (time 
lagged) approach. This explicit approach was proved to be sufficient for shape change coupling simulations 
demonstrated in Refs. 1 and 2. However, for explicit coupling simulations with surface thermal chemistry included, 
numerical oscillation of predicted surface temperature was observed at conditions with relatively high surface heat 
fluxes. 5 The time step size for updating the interface boundary conditions had to be reduced to minimize the 
numerical oscillation. Consequently, the total computational time required to complete a simulation increased 
significantly. Additionally, for some explicit coupling simulations, the numerical oscillation of surface temperature 
became so severe that a converged solution was not obtained. Thus, an alternative approach needs to be developed 
to efficiently perform fluid-solid surface thermal chemistry coupling at high surface heating environments. 

The purpose of this paper is to demonstrate an implicit DPLR/3dFIAT coupling simulation system and its 
applications to the arc-jet simulation and the trajectory-based simulation of Stardust Earth Reentry Capsule. The 
boundary condition with general non-equilibrium finite-rate chemistry for gas/surface interactions implemented in 
the DPLR code by MacLean 6 is used to predict char mass injection rate. In this simulation system, the surface 
species mass balance is performed in DPLR, and the surface energy balance is performed in 3dFIAT. The chemical 
equilibrium assumption typically used in material thermal response code can be removed. The hot-wall ablating 
convective heat flux is directly computed in DPLR based on the surface temperature and pyrolysis gas injection rate 
computed in 3dFIAT. The non-equilibrium gas/surface interaction chemistry between air and carbon surface is 
based on that developed by Park 7 and is modified to better match the arc-jet data. Coupled fluid-material response 
analyses of stagnation tests conducted in NASA Ames Research Center arc-jet facilities are first performed for code 
validation. The ablating material used in these arc-jet tests is a Phenolic Impregnated Carbon Ablator (PICA). 8 The 
predictions of using implicit coupling algorithm are compared with arc-jet data for seven selected CEV ADP arc-jet 
tests. Then, the implicit approach is used for trajectory-based simulation. A trajectory based simulation for Stardust 
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Earth Reentry Capsule is computed successfully using the implicit approach. The predictions are presented and 
discussed in detail. 


II. Surface Thermal Chemistry and Shape Change Coupling 

The code coupling is basically the management of information exchange to satisfy all the conservation laws at 
the interfaces. Shape change coupling tracks the time-dependent moving interfaces, and thermochemistry coupling 
enforces the mass and energy balances at the interfaces. Shape change coupling also requires a moving volumetric 
grid system to reflect the corresponding moving interfaces. Thermochemistry coupling updates the boundary 
conditions on both sides of interface. If the boundary conditions are updated in a time-lagged manner, it is called an 
explicit approach. If the boundary conditions are updated as part of the solution, then it is called implicit approach. 
For the implicit approach in this paper, the interface boundary conditions are assumed to be piecewise linear in time, 
and additional iterations are performed to obtain the updated interface boundary conditions. 

The interface boundary conditions between a TPS material and its surrounding flow field can be defined by 
solving species mass conservation and energy balance equations. Species mass conservation at the surface of TPS 
material is written as: 9 

-pD i VC i + pv w C t = M- w z + m g C ig (1) 

The first term on the left-hand side is mass transfer through diffusion, and the second term is mass transfer due to 
convection. On the right-hand side are the source terms due to gas-surface interaction and pyrolysis gas injection. 
Based on global mass balance at the surface, the following equation for the total mass blowing rate is expressed as: 

pv w =m c +m g (2) 

The total convective heat flux to the surface for the flow field that includes a two-temperature model is given as: 

Vconv = -k t VT t - k v VT v + Y J j i,P D i S JC i . (3) 

Energy conservation equation at the surface of TPS material is written as 

Qconv "f "f TTlgijXg T OL w C[ rw <J£ W (T W Too) Qcond ^ (4) 

The first term in Equation (4) is the total convective heat flux, the second and third terms represent the heat of 
ablation, the fourth and fifth terms are radiation absorption and emission, respectively, and the final term is the rate 
of heat conduction into the TPS material. 

To obtain the solutions of Equations (1) to (4) requires the computations of thermal and species diffusion rates of 
flow field at the surface as well as thermal diffusion and pyrolysis gas injection rates of TPS material at the surface. 
The governing equations for both fluid and solid can be solved simultaneously along with Equations (1) through (4). 
However, for many simulations, flow solver and material response code are two independent programs, and a 
coupled simulation is usually applied. In a coupled simulation, the governing equations for fluid and solid are solved 
separately, and Equations (1) to (4) are solved either in the flow code or in the material thermal response code. Thus, 
communication between two codes needs to be established for exchanging information on surface thermal chemistry 
and shape change. 

If the material surface is at chemical equilibrium, both the mass and energy conservation equations at the surface 
can be performed in the material simulation code as described in our previous work. 1,2 For a chemical equilibrium 
surface, the chemical species at the surface are determined by using a chemical equilibrium code, such as ACE 10 or 
MAT. 11 In our chemical equilibrium analyses, the hot wall ablating heat flux was estimated from the cold wall non- 
ablating heat flux using an engineering correlation with a blowing reduction parameter. Also, species mass transfer 
rate was assumed to be proportional to heat transfer rate based on a constant Lewis number. Under these 
assumptions, surface thermal chemistry coupling is not required, and only shape change information has to be 
exchanged between flow solver and material response code to correctly predict the aerothermal environments and 
material thermal response. 

For a general finite-rate surface boundary condition, the approach taken here is solving the species mass 
conservation equation (Eq. 1) with the flow- field governing equations and solving the total energy balance equation 
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(Eq. 4) with the solid material governing equations. If the information of surface thermochemistry and shape change 
shared between flow solver and material response solver is time lagged, it is an explicit approach. For explicit 
coupling approach, the interface boundary conditions for material response code are assumed to be piecewise 
constant in time. In other words, the boundary conditions during At D are equal to those at t*. The time step to update 
the interface boundary conditions, At D , can be varied dynamically. This time step may become unrealistically small 
due to numerical instabilities. For implicit coupling approach studied in this paper, the boundary conditions for 
material thermal response code are assumed to be piecewise linear in time. The boundary conditions during each At D 
are calculated by linear interpolation between t* and t*+ At D . The interface boundary conditions at t*+ At D are 
unknown and thus have to be obtained through iterations. This approach is more computationally time intensive as 
compared with explicit coupling for completing each time cycle of At D . However, the implicit approach can reduce 
the numerical oscillations, and may decrease the total run time. If the surface boundary condition is a stiff function 
of time or surface temperature, such as that for a sublimating carbon surface, the implicit approach is superior to 
explicit approach for conducting an integrated fluid-solid simulation. Either a constant time step or a variable time 
step can be used for updating the interface boundary conditions implicitly, as long as time accuracy is maintained. 

The schematic diagram in Figure 1 depicts how the implicitly coupled simulation is performed. This is a time 
accurate simulation, and computations start from a time marching in the material response simulation code using the 
boundary conditions of cold-wall non-ablating convective heat flux and surface pressure estimated based on initial 
free stream conditions. Time integration is carried from t = t* to t = t*+ At D . Since the boundary conditions at the 
end of integration are unknown, they have to be obtained through iterations. In the first iteration, the boundary 
conditions at t = t*+ At D is set to be equal to those at t = t*. This is the initial guess of the interface boundary 
conditions at t*+ At D . At the end of each time integration, if the maximum surface temperature variation at t = t*+ 
At D exceeds the pre-determined value, the material thermal response computation is temporarily put in a waiting 
mode and flow field simulation is performed using the latest predicted surface temperature and pyrolysis gas 
injection. The pyrolysis gas is assumed to be at chemical equilibrium before being injected into the adjacent air. This 
assumption was proved to be reasonable for conditions studied in this paper. 12 The computational grid system for the 
flow field is reconstructed based on the shape change predicted by the material response code. The free stream 
conditions are also updated according to the proposed flight trajectory. Each flow simulation is a steady-state 
computation. When the steady-state flow/radiation solution is obtained, the time-dependent material response 
simulation is resumed using the newly predicted surface heat flux, pressure, char mass injection rate, and surface 
radiation as the boundary conditions at the time step of t*+ At D . It usually takes about 3 to 5 cycles of iterative loop 
to obtain the converged boundary conditions at t*+ At D . If there is no iterations performed, this approach converts to 
the explicit coupling. Thus, the explicit approach is just a simplified version of the implicit approach. The explicit 
approach is similar to the implicit approach in which only one single cycle of iterative loop is performed during each 
time marching interval of At D . This iterative process repeats for each At D until the end of flight trajectory or arc-jet 
exposure. 
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Figure 1 Flow chart for implicit surface thermal chemistry and shape change coupling. 


III. Chemistry Models 

Two sets of non-equilibrium gas/surface interaction models between air and a carbon surface were found in the 
open literature. The first set was developed by Park 7 . In our previous study, the nitridation reaction in the original 
Park’s model was replaced by the nitrogen recombination reaction to be consistent with observations in arc -jet tests 
using pure nitrogen gas. 5 The second set was developed by Zhluktov et. at. 13 However, it was found in our previous 
work 5 that Zhluktov’ s model far under-predicted PICA recession. Thus, Zhluktov’ s model is not considered in this 
work, and Park’s carbon sublimation model is revised to better match the arc -jet data. The modified Park’s 
gas/surface reactions used in this paper are listed below: 

Modified Park Carbon/ Air Chemistry: 

#P1 O + C(b) -> CO 
#P2 0 2 + 2C(b) -> 2CO 
#P3 2N -> N 2 
#P4 3C(b) -> C 3 
#P5 C 3 -> 3C(b) 

Here C(b) represents the solid carbon species. In reaction #P3, nitrogen atoms are assumed to be fully 
recombined at the surface. To better agree with arc-jet data in the sublimation region, the pre-exponential factors of 
reaction rate constants for reactions #P4 and #P5 used in this work are 8.27x1 0 14 and 0.25. The rest of rate constants 
can be found in Refs 6 and 9. Marschall developed a general formulation for finite-rate gas and surface 
interactions. 14 The DPLR code was enhanced by MacLean based on Marschall’ s work to solve the surface species 
mass balance equation for the finite-rate gas surface reactions. The detail of his implementation can be found in Ref. 
6. All flow simulations presented in this work are performed using this version of DPLR code (DPLR/4.03.0cvs). 

There are twenty nine gas phase chemical species used in this study for the simulation of PICA and air- 
Argon interactions. The chemical species are C0 2 , CO, N 2 , 0 2 , NO, C 2 , C 3 , CN, H 2 , HCN, C, N, O, H, CH, CH 2 , 
C 2 H, C 2 H 2 , C 4 , C 3 H, C 4 H, Ar, N2 + 02 + NO + C + N + 0 + and e". These species were selected based on a chemical 
equilibrium computation for a PICA/air-Argon system. Their enthalpy makes up more than 95% of total gas 
enthalpy compared with the baseline 119 chemical equilibrium species model developed by Orion TPS Advanced 
Development Project. 15 The gas phase chemical reactions considered in the simulations are as follows. 

#1 co 2 + m~co + o + m 

#2 CO + M<->C + 0 + M 
#3 N 2 + M~N + N + M 
#4 0 2 + M~0 + 0 + M 

#5 NO + M<->N + 0 + M 
#6 C 2 + M~C + C + M 
#7 C 3 + M<->C 2 + C + M 
#8 CN + M<->C + N + M 
#9 H 2 + M~H + H + M 
#10 N 2 + Ar<->N + N + Ar 
#11 0 2 + Ar <-> O + O + Ar 
#12 NO + Ar <-> N + O + Ar 
#13 N0 + 0~0 2 + N 
#14 N 2 + 0~N0 + N 
#15 C0 + 0~0 2 + C 
#16 C0 2 + 0~0 2 + C0 
#17 C0 + C~C 2 + 0 
#18 CO + N <-> CN + O 
#19 N 2 + C~CN + N 
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#20 CN + O <-> NO + C 
#21 CN + C~C 2 + N 
#22 HCN + H~CN + H 2 
#23 CH + M<->C + H + M 
#24 CH 2 + M~C + H 2 + M 
#25 CH 2 + M~CH + H + M 
#26 C 2 H + M~C 2 + H + M 
#27 C 2 H 2 + M~C 2 H + H + M 
#28 C 2 + C 2 ~C 3 + C 
#29 C 2 + H 2 ~C 2 H + H 
#30 CH 2 + C~CH + CH 
#31 CH 2 + CH 2 ^ C 2 H 2 + H 2 
#32 CH 2 + C~C 2 H + H 
#33 CH 2 + C 2 H <-> CH + C 2 H 2 
#34 CH + CH~C 2 H + H 
#35 CH + C 2 H <-> C 2 H 2 + C 
#36 CH 2 + CH 2 ~C 2 H 2 + H + H 
#37 C 2 + C 2 + M~C 4 + M 
#38 C + CH~C 2 + H 
#39 C + C 2 H~C 3 + H 
#40 C 2 + CH~C 3 + H 
#41 C 2 + C 2 H~C 4 + H 
#42 CH + H~C + H 2 
#43 CH 2 + H~CH + H 2 
#44 C 2 H + C 2 H <-> C 2 H 2 + C 2 
#45 C + C 2 H 2 ~C 3 H + H 
#46 C 2 + C 2 H 2 <-> C 4 H + H 
#47 C 2 + C 4 H ~ C 2 H + C 4 
#48 C 2 H + C 2 H <-> C 4 H + H 
#49 C 3 H + H~C 3 + H 2 
#50 C 4 H + H~C 4 + H 2 
#51 H + C 4 + M~C 4 H + M 
#52 CH + CH~C 2 + H + H 
#53 C 2 H + H 2 <-> C 2 H 2 + H 
#54 C + e <-> C + + e + e 
#55 N + e" <-> N + + e" + e" 

#56 O + e <-> 0 + + e + e 
#57 N + O <-> NO + + e" 

#58 N + N <-> N 2 + + e" 

#59 O + O ~ 0 2 + + e 
#60 N + + N 2 ~N 2 + + N 
#61 0 + + N 2 ~N 2 + + 0 
#62 0 2 + + O <-> 0 + + 0 2 
#63 0 + + N0~N + + 0 2 
#64 NO + + 0 2 <-> 0 2 + + NO 
#65 NO + + N ~ N 2 + + O 
#66 NO + + O ~N + + 0 2 
#67 0 2 + + N <-> N + + 0 2 
#68 0 2 + + N 2 ^ N 2 + + 0 2 
#69 NO + + N ^ 0 + + N 2 
#70 NO + + 0 <-> 0 2 + + N 

The reaction rates for reactions #1-22 are taken from the work of Olynick et al. 16 for Stardust earth entry 
simulations. The rates for reactions #23-36 are taken from Gokgen’s 17 paper for simulations of Titan atmosphere 
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entry. The rates for reactions #37-53 are taken from the study of Kruse et al. 18 for high-temperature pyrolysis of 
Acetylene. The rates for reactions #54-70 are taken from Park et al. 19 for Mars entries. 


IV. Results 

In this section, we present two sets of computation performed using the DPLR/3dFIAT coupled simulation 
system discussed in the previous sections. In the first set of computation, arc-jet simulations are performed. This is 
to demonstrate that the implementation of implicit coupling technique is self consistent and the time accurate 
solutions can be achieved. Predictions are compared with those predicted using the explicit approach, and with the 
data from seven selected Crew Exploration Vehicle Advanced Development Program (CEV ADP) arc-jet tests. In 
the second set of computations, trajectory-based simulations are performed for Stardust fore-body heat-shield during 
Earth reentry. Attempt has been made using the explicit approach to perform the same trajectory based simulation, 
but computation failed to converge around the peak heating point due to excessive numerical oscillation in predicted 
surface temperature. 

A. Arc-Jet Simulation 

Coupled DPLR/3dFIAT computations using the revised Park’s finite-rate air-carbon gas/surface interaction 
model described early in this paper are performed for arc-jet flow over a stagnation test model. Model validation is 
accomplished by comparing the predictions with the data from seven selected CEV ADP arc-jet tests conducted over 
a range of stagnation heat flux and pressures ranging from 107 W/cm 2 at 2.3 kPa to 1100 W/cm 2 at 84 kPa. 20 The 
test conditions for these seven test cases are listed in Table 1. 


Table 1: Seven arc-jet cases selected for detailed analysis. 


Case number 

Stagnation point heat flux, 
W/cm 2 (cold wall) 

Stagnation point pressure, 
kPa 

Exposure time, 
s 

Total enthalpy 
MJ/kg 

1 

107 

2.3 

55 

15.2 

2 

169 

5.0 

60 

17.0 

3 

246 

8.5 

42 

19.3 

4 

395 

17.2 

34 

21.4 

5 

552 

27.3 

30 

23.3 

6 

694 

31.0 

27 

29.2 

7 

1102 

84.4 

10 

25.6 


Figure 2 shows the geometry and material map of stagnation model used in these tests. The model nose radius is 
equal to the model diameter (r n = r D = 5.08 cm, and r c /r n = 1/16), and the sides are cylindrical. This is the so called 
“iso-q” geometry. The tested TPS material is PICA. The detailed comparisons between two finite-rate air-carbon 
gas/surface interaction models (Park’s and Zhulktov’s) using a time lagged (explicit) coupling were discussed in our 
previous paper. 5 In this study, the simulations for the same seven arc -jet cases are performed using the implicit 
approach. The PICA material in-depth thermal response and flow field structure predicted using implicit coupling 
are not presented in this paper, because they are essentially the same as those predicted using the explicit approach, 
and were discussed extensively in our previous paper. 2 Here, we focus on the issue of numerical oscillation 
associated with fluid/solid coupled simulations, and the advantage/disadvantages of applying the implicit coupling 
to update the interface boundary conditions. Three of these seven arc-jet cases will be discussed in detail to explore 
the difference between the predictions using implicit and explicit approaches. Case 2 represents the low end side of 
heat flux (169 W/cm 2 ), while Cases 6 and 7 represent those at the high end side (694 and 1 102 W/cm 2 ). 
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PICA Air Gap LI-2200 


Figure 2 Geometry and material map of the model used in the arc-jet tests. 

Figure 3a presents the predicted stagnation point temperature history for Case 2 during the first 5 seconds of 
simulations using both implicit (black line) and explicit (red line) coupling approach. Based on our experience, 
numerical oscillations in the surface temperature and mass injection rate normally occur during the first few seconds 
of arc-jet exposure because of the abrupt change to the surface environment. Thus, the comparison between these 
two predictions after the first 5 seconds is not presented. Since the surface heat flux is relatively benign for this case, 
the maximum stagnation point temperature is around 2200 K. The difference on predicted surface temperature 
between two coupling approaches is small enough to be ignored. Figure 3b shows the comparison of stagnation 
point mass injection rates of char and pyrolysis gas. The pyrolysis gas mass injection rate is higher than the char 
mass injection rate during the first 5 seconds. Aside from a few minor oscillations in the first few time steps of 
explicit simulation, the agreement is considered to be very good between these two approaches. In the implicit 
simulation, a fixed time step, At D , of 0.1 s is used. This is the time step for DPLR and 3dFIAT to exchange interface 
information, not the time marching step in 3dFIAT. For the explicit approach, a variable time step is adopted to 
minimize the total computational time. In the first few time steps, the step size is of the order of 10" 2 s, and then 
gradually it increases up to about 0.25 s. The time step size used in these two coupling simulations are quite 
different, but the predicted interface conditions, including temperature and mass injection rate, are almost identical. 
Since the numerical oscillations associated with the explicit coupling is insignificant under this surface heating 
condition, the explicit coupling with variable time step size appears to be more time efficient. The computational 
time for each implicit time step is equivalent to 3 to 5 explicit time steps, due to additional iterations. Unless the 
time step for the implicit simulation is 3 to 5 times greater than that of the explicit simulation and the time accuracy 
can still be maintained, the explicit coupling should be adequate. Thus there is no advantage of applying the implicit 
coupling for this case. 

Similar comparison for Case 6 is shown in Figs 4. Case 6 has stagnation point heat flux of 694 W/cm 2 , which is 
much higher than that of Case 2. Figure 4a is the comparison of stagnation point temperature history and Fig. 4b 
presents stagnation point mass injection rate history during the first 5 seconds of arc-jet exposure. The maximum 
surface temperature reaches around 3100 K. The variable time step is adopted in both computations. Numerical 
oscillations are observed in the results predicted by both approaches. The predictions using implicit approach (black 
lines) have much less variation in amplitude than those using explicit approach (red lines). Regardless of numerical 
oscillations, the agreement on predicted interface boundary conditions between these two approaches is reasonably 
good. During the first second of computation, the time step size used in the explicit coupling (order of 1 0" 3 s) is far 
smaller than that of the implicit coupling (order of 10" 2 s), but numerical oscillations still cannot be further reduced. 
The predictions after the first 5 seconds have no numerical oscillation and thus are not shown in Figs. 4. To further 
reduce numerical oscillations in the predictions using the explicit approach is impractical. Since the required time 
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step size becomes unreasonably small. The implicit approach thus is a better choice for simulations for cases similar 
to that of Case 6. The total run time of Case 6 for the implicit coupling is about 25% of that of the explicit coupling. 
The similar comparisons between the implicit and explicit approaches for Case 7 are shown in Fig. 5. The time step 
used for the explicit coupling is about one order of magnitude less than that of Case 6, while that for the implicit 
coupling remains the same as Case 6. The numerical oscillations are similar to those observed in Case 6. For 
conditions with even higher heat fluxes, the explicit coupling usually fails to converge because of excessive 
numerical oscillations even if the selected variable time step size is further reduced. The implicit approach thus 
becomes the more practical way to perform simulations at such conditions. Generally speaking, the agreement on 
predicted interface boundary conditions between implicit and explicit simulations is good as long as the numerical 
oscillation is moderate. The implicit approach has no instability issues but may have minor numerical oscillations, 
which mostly can be minimized by controlling the time step size. A fixed constant time step can be used in the 
implicit simulations for most applications. A variable time step is only used to minimize the numerical oscillations. 
The explicit approach requires less computational time to complete one single time step for updating the interface 
conditions, and its coding is fairly straightforward as compared with the implicit approach. The best application for 
the explicit approach is the fluid-solid chemical equilibrium shape change coupling, in which surface temperature 
and surface thermal chemistry are not involved in the coupling process. 

The predicted stagnation point total recession values for all seven arc-jet cases are presented in Fig. 6, and are 
compared with the measurements. As mentioned in the previous section, the pre-exponential factors used in our 
surface sublimation reactions, #P4 and #P5, are different from those suggested by Park. The pre-exponential factors 
suggested by Park under-predicted the carbon sublimation by about 20%. 5 The maximum difference between the 
computed and measured recession values for all seven cases is less than 5 %. The comparison of predicted 
stagnation point ablating hot-wall heat flux for the non-equilibrium and equilibrium surface chemistry model is 
shown in Fig. 7. The non-equilibrium model is the current finite-rate gas/surface interaction model, where as the 
chemical equilibrium model is our base-lined model with blowing reduction parameter of 0.5. The predictions using 
the chemical equilibrium surface assumption are just slightly lower than those using the finite-rate model for all 
seven cases. Thus, this indicates that the chemical equilibrium surface is a reasonable approximation for prediction 
of stagnation point heat flux and surface recession under arc stream conditions studied in this section. 
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Figure 3 Comparison of stagnation-point time histories of various surface quantities for Case 2 using the 
implicit and explicit coupling approaches. 




Figure 4 Comparison of stagnation-point time histories of various surface quantities for Case 6 using the 
implicit and explicit coupling approaches. 
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Figure 5 Comparison of stagnation-point time histories of various surface quantities for Case 7 using the 
implicit and explicit coupling approaches. 
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Figure 6 Comparison of the computed stagnation-point 
total recession with the arc-jet data. 


Figure 7 Comparison of the computed stagnation-point 
hot-wall heat fluxes for the chemical non-equilibrium 
and equilibrium models. 
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B. Trajectory Based Simulation - Stardust Earth Reentry Capsule 

The predictions of trajectory based simulation using implicit coupling between DPLR and 3dFIAT for Stardust 
fore -body heat-shield is presented in this section. The Stardust fore-body heat-shield material is PICA. The PICA 
surface was expected to reach the low sublimation region around the peak convective heating point. Thus, this is an 
ideal case to demonstrate the implicit coupling approach. The surface geometry for Stardust fore-body heat shield is 
an axi-symmetric 60-degree half angle sphere-cone with a nose radius of approximately 23 cm. The size of 
computational grids for flow field calculation is 81 x 91, and that of material thermal response simulation is 81x61. 
The finite-rate gas phase chemistry and gas/surface interactions used in this simulation are the same as those used in 
the arc-jet simulation discussed in the previous section. In arc-jet simulation, the arc stream conditions are fixed in 
time, and the radiation emitted from the arc stream is negligibly small. For the trajectory-based simulation discussed 
in this section, the NEQAIR2009 code is used to estimate the radiation emitted from shock layer to the surface. 21 
The trajectory used in this simulation is shown in Fig. 8. 22 An attempt has been made using the explicit approach to 
perform this simulation. However, because of sever numerical oscillations starting at around t = 50 s, no converged 
solution was obtained. Thus, comparison of predictions using the implicit and explicit coupling approaches is not 
available. 

Similar trajectory based simulation was performed in Ref. 16, in which the gas/surface coupling simulation was 
conducted at 11 pre-selected trajectory points from 34 to 80 s. Chemical equilibrium assumption was made at the 
heat-shield surface, and the effect of shape change was neglected. Additionally, the one-dimensional FIAT 23 code 
was used for material ablation and thermal response simulation, and the NOVAR code was used to predict shock 
layer radiation in the work of Olynick et al. 24 

Before t = 35 s, the solutions of full Navier-Stokes equations are not computed, due to low density effect. The 
continuous flow assumption may be in question. After 100 s, the obtained Navier-Stokes solutions from DPLR were 
not fully converged. The surface convective heat flux at approximately 20 s is assumed to be equal to 1 W/cm 2 , and 
that between 20 and 35 s is estimated by linear interpolation. 3dFIAT starts the material response simulation at t = 0 
s, and is not coupled with DPLR and NEQAIR to exchange interface boundary conditions until 35 s. Implicit 
coupling approach is used starting from 35 s and ending at 100 sec. The time step size for implicit coupling, At D , is 
set to a constant of 1 s. In other words, every one second the interface boundary conditions between gas and carbon 
surface have to be updated implicitly. It usually requires 3 to 5 internal iterations to converge. 

The predicted stagnation point heating history is shown in Fig. 9a. The maximum hot wall ablating convective 
heat flux is about 600 W/cm 2 and the maximum radiation is about 90 W/cm 2 . Both convective and radiative heat 
fluxes at the stagnation point predicted by Olynick (symbols) are also presented in the same chart for comparison. 
Olynick’ s convective heat flux is about 11% higher than that predicted by current DPLR/3dFIAT simulation system. 
The bifurcation diffusion model and chemical equilibrium surface assumption adopted in Olynick’ s work tend to 
predict higher mass diffusion rate, and surface heat flux. The non-ablating cold wall heat flux (dashed line) is also 
shown in this chart. The maximum convective cold wall heat flux is slightly above 900 W/cm 2 . The stagnation point 
mass injection rates (solid lines) and surface temperature (dashed line) histories are presented in Fig. 9b. The 
maximum char mass injection rate is slightly above 0.05 kg/m 2 -s. The predicted maximum char mass injection rate 
in Olynick’s work is approximately 0.08 kg/m 2 -s. The pyrolysis gas injection rate remains fairly uniform for a long 
period during the trajectory, and the maximum value is near 0.01 kg/m 2 -s. The predicted maximum surface 
temperature reaches about 3300 K at the peak convective heating point. 

Figures 10a and 10b show the convective and radiative heat fluxes along the fore-body surface at t = 40, 50, 56, 
62, and 72 s. The heating distributions are similar for all these time points. The maximum heating is located at the 
stagnation point, a sharp decline on surface heating is observed across the junction between sphere and cone, and a 
spike occurs along the surface of heat-shield corner. The peak radiative heating occurs at around t = 50 s, which is 
not at the same time as that of peak convective heating. 

Char and pyrolysis mass injection rates along fore-body surface at time points equal to 56 and 72 s are shown in 
Figs. 11a and 1 lb, respectively. The distribution of mass injection rate of pyrolysis gas shown in Fig. 11a appears to 
be fairly uniform over most of the sphere-cone surface except at the corner region. The char mass injection rate 
declines rapidly along the spherical section, and then declines moderately in the conical section, as shown in Fig. 
lib. 

Figures 12a and 12b present species mass fractions along stagnation stream line at two time points, 56 and 72 s. 
At 56 s, the free stream enthalpy is high enough to fully dissociate nitrogen molecules within the shock layer. At 76 
s, the dissociation rate of nitrogen molecules is very low because of much lower total free stream enthalpy compared 
with that of 56 s. The major ablation product is CO at both time points. At 56 s, the carbon surface is in the low 
sublimation region. Thus, the second most abundant ablation product is C 3 . The third one is hydrogen atom and the 
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forth one is hydrogen molecule. At 72 s, the second most abundant ablation product is H 2 , which is coming from 
pyrolysis gas, and the H mass fraction is one order of magnitude lower than H 2 . The C 3 mass fraction is very low, 
because the carbon surface reaction is in the regime of diffusion controlled oxidation. 

Figure 13 presents the stagnation point recession history. The predicted total recession is approximately 7 mm at 
100 s, which is about 33% less than that predicted using the chemical equilibrium surface assumption described in 
Olynick’s paper. The post flight measurement of surface recession at the stagnation point is around 6 mm. The 
stagnation point recession is still over-predicted by about 16% by using current finite-rate coupled simulation 
system. 
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Figure 8 Flight trajectory for the Stardust Earth 
reentry capsule. 





b) Surface temperature and pyrolysis and char mass 
a) Ablating hot-wall convective and radiative heat injection rates, 
fluxes. 


Figure 9 Computed time history of the stagnation point surface quantities. 
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a) Convection b) Radiation 

Figure 10 Heat flux distributions along fore-body heat-shield at various times on the trajectory 




a) t = 56 s b) t 72 s 


Figure 11 Mass injection rate distributions along the fore-body heat-shield at two trajectory points. 
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a) t = 56 s 


b) 72 s 


Figure 12 Species mass fractions along the stagnation streamline at two trajectory points 



Figure 13 The computed stagnation point total recession 
as function of time. 
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V. Conclusions 

In this study, coupled simulations of a material thermal response code and a flow solver with non-equilibrium 
gas/surface interaction for simulation of charring carbon ablators are performed using either implicit or explicit 
approaches. For analyses of the stagnation tests conducted in NASA Ames Research Center arc-jet facilities both 
implicit and explicit approaches are used. It is found that the implicit approach is a better approach at the high 
heating conditions, while the explicit coupling is more time efficient at the lower heating conditions. The implicit 
approach requires more computational time to complete a single time step for updating the interface boundary 
conditions due to additional iterations. However, the implicit coupling appears to have no numerical instability issue 
and does not have any severe numerical oscillations for all seven test conditions. On the other hand, the explicit 
approach is simpler to construct, but it may require very long run time or fail to converge if excessive numerical 
oscillations in the predicted surface temperature do occur. The predicted surface recession using the modified Park’s 
model was compared with the experimental measurements at stagnation cold- wall heat fluxes ranging from 107 to 
1100 W/cm 2 . The maximum difference between predicted recession and arc-jet data is less than 5%. 

Trajectory-based simulations for Stardust fore-body heat-shield was also performed using the implicit approach. 
A constant time step is adopted for this case. The computations using the explicit coupling were not successful 
because of numerical oscillations around the peak heating point. It is shown that the implicit approach is more 
appropriate for a trajectory-based surface thermochemistry coupling for Stardust entry. Using the current best 
available surface chemistry, the predicted stagnation-point heat flux and recession for the Stardust capsule are 
smaller than those predicted using the chemical equilibrium surface assumption by 11% and 33%, respectively. 
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